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We propose a method to probe the nature of phase transitions in lattice QCD at finite temperature 
and density, which is based on the investigation of an efi'ective potential as a function of the average 
plaquette. We analyze data obtained in a simulation of two-flavor QCD using p4-improved staggered 
quarks with bare quark mass m/T = 0.4, and find that a first order phase transition line appears in 
the high density regime for Hq/T ^ 2.5. We also discuss the difference between the phase structures 
of QCD with non-zero quark chemical potential and non-zero isospin chemical potential. 



. . . PACS numbers: 11.15.Ha, 12.38.Gc, 12.38.Mh 

o 

O . I. INTRODUCTION 

CN 

^ ' In the last several years remarkable progress has been made in numerical studies of lattice QCD at finite temperature 
/"N , (T) and quark chemical potential (/iq). The transition line, separating hadron phase and quark-gluon plasma (QGP) 
" ' phase, was investigated from /x, = to finite Hq P, H, S, Ij, [![ , and the equation of state was also computed at low 
density [1,11,0,11,11. Among others, the study of the endpoint of the first order phase transition line in the (T, fiq) 
plane is particularly important both from the experimental and theoretical point of view. This existence of such a 
critical point is suggested by phenomenological studies [H, [ll|, [ll| . The appearance of the critical endpoint in the 
(T, fiq) plane is closely related to hadronic fluctuations in heavy ion collisions and may be experimentally examined 
by an event-by-event analysis of heavy ion collisions. 

Although many trials have been made to prove the existence of the critical endpoint by first principle calculation 
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^ ' in lattice QCD, no definite conclusion on this issue is obtained so far. The first trial to find the critical endpoint by 
, numerical simulations was performed in Ref. investigating the finite size scaling behavior of Lee- Yang zeros in the 
complex /? = 6/(7^ plane. The difficulty in the Lee- Yang zero method for finite density QCD is discussed in Ref. [11]. 
The radius of convergence in the framework of the Taylor expansion of the grand canonical potential can establish a 
^ . lower bound on the location of the critical endpoint [3, [31 • There are also studies in which the behavior of the 
^\ ' critical endpoint as a function of the quark masses is examined by using the property that a critical endpoint exists at 
■rj" . fiq = in the very small quark mass region for QCD with three flavors having degenerate quark masses fisl [l6l. Il7t . 
■ Moreover, studies by simulations of phase-quenched finite density QCD, have been performed in Ref. H, ig. l20j|. 
" The purpose of this study is to clarify the existence of the endpoint of the first order phase transition line in 
\^ ■ the {T,ijq) plane. We propose a new method to investigate the nature of transition. In the study of finite density 
' lattice QCD, the reweighting method [2l|, [2^ plays an important role. However, the calculation of physical quantities 
C . becomes increasingly more difficult for large fiq due to the sign problem [H, [2^ . We also consider a way to avoid the 
sign problem. 

We evaluate an effective potential as a function of the average plaquette, and identify the type of transition from 
the shape of the potential. The partition function can be written as 



Z{p,fiq) = J R{P,^,q)w{P)e-''^^P'^^ dP, (1) 

where P denotes the plaquette value, Sg{P,(3) is the gauge action, w{P) is the state density at fiq ~ for each P, 
and R{P, fiq) is the modification factor for finite fiq. R{P, fiq) is obtained by calculating the quark determinant det M 
and is assumed to be real and positive. We then define the effective potential as V{P, (3, fiq) — — \n{Rwe~^!>). If there 
is a first order phase transition point, where two different states coexist, the potential must have two minima at two 
different values of P. However, the calculation of the quark determinant is quite expensive and is actually difficult 
except on small lattices. Moreover, the sign problem is serious when we calculate R{P^ fiq) directly. 

This study is based on the following two ideas to avoid these problems. One is that we perform a Taylor expansion 
of Indet M (fiq) in terms of fiq around fiq = and calculate the expansion coefficients, as proposed in Ref. [3|. The 
Taylor expansion coefhcients are rather easy to calculate by using the stochastic noise method. Although we must cut 
off this expansion at an appropriate order in fiq, we can estimate the application range where the approximation is 
valid for each analysis. While the application range of the Taylor expansion of InZ should be limited by the critical 
point because InZ is singular at the critical point, there is no such limit for the application range in the expansion of 
In i?(P, fiq) because the weight factor should always be well-defined. This discussion is given in Sec. IIIIBI 
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The second idea is that we consider the probabiHty distribution function in terms of the complex phase of the quark 
determinant 9 when P and | det M\ are fixed. We assume the distribution function is weU-approximated by a Gaussian 
function, and perform the integration over 0. If we adopt this assumption, the sign problem in the calculation of 
In i?(P,/ig) is completely solved. This assumption is reasonable for sufficiently large volume and is suggested by the 
simulation results given in this study. We discuss this method in Sec. IIII CI 

General remarks on the phase transition in lattice QCD are given in Sec. |TT1 and an effective potential as a function 
of the average plaquette is introduced. We discuss the reweighting method for the study of the QCD phase structure 
at non-zero temperature and density in Sec. IIIII We evaluate the effective potential using data obtained with two- 
flavors of p4- improved staggered quarks in Ref. [7|. We also discuss the phase structure of QCD with isospin chemical 
potential. Our conclusions are given in Sec. IIVI 



II. PROBABILITY DISTRIBUTION FUNCTION AND PHASE TRANSITION 



The grand canonical partition function of lattice QCD is given by 

Z{(3,fig)^ I VU {det Mf'e~^^, (2) 

and the expectation value of an operator O is calculated by 

(O) = ^ J VUO{detM)^'e-^\ (3) 

where M{^q) is the quark matrix. Ni is the number of flavors. When one use a staggered type quark action, TVf is 
replaced by Nf/4. Sg{l3) is the gauge action, which is given by a linear combination of the Wilson loops W^^"' (x), 
where I y. J , jiv and x are the size, direction and position of the Wilson loop, respectively, /3 is a simulation parameter 
related to the gauge coupling g being j3 = &/g^- The simplest gauge action is the standard plaquette action given by 
the following equation, 

s, = -pY1 K^'i^)- (4) 

Because the 1x1 Wilson loop is defined on an elementary square (plaquette), W^^^ is usually called plaquette or 
plaquette variable. 

In a Monte Carlo simulation, we generate configurations of link variables {J7^(a;)} with the probability in proportion 
to the weight factor (det M)^' e"'^^ and the state density of {[/^(x)}. The expectation value is then estimated by 
taking an average of the operator 0[U^] over the generated configurations {U^{x)}. 

We introduce a probability distribution function of the plaquette, w{P), which is defined by 

w{P') = Jt^U 5{P' - P) (det M)^f e^^^^-'^^, (6) 

where 5{x) is the delta function. For later discussions, we define the average plaquette P as P = —Sg/{6f3Nsitc)- This 
is the average of the plaquette over all elementary squares for the standard gauge action, Eq. ^ . NsHc = x Nt is 
the number of sites. Using the distribution function, the expectation value can be rewritten as 

{0[P])ip) = ^J 0[P] w{P) dP, Z = j w{P) dP, (7) 

for an operator given by the plaquette ^^[P]. In the calculation of Eq. ([6|), we actually use an approximate delta 
function such as a box type function, S{x) « (for A/2 < x < A/2),0 (otherwise)}, or a Gaussian function, 

S{x) « l/(A-y/7r) exp[— (x/A)^]. For the case of the box type, we can estimate w{P) by counting the number of 
configurations for each value of P with the width of box A. As A decreases, the approximation becomes better but 
the statistical error becomes large because the number of configurations in each block becomes small. Hence, we must 
adjust the size of A appropriately. 
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Next, we discuss the shape of the probabihty distribution function. In general, the number of states increases 
exponentiaUy as the gauge fields become random. On the other hand, the random configurations are exponentially 
suppressed by the weight factor exp(6/3-/Vsitc-P), since the plaquette is one when the gauge field is Ufj,{x) = 1 uniformly 
(free gas limit) and P decreases as the configuration becomes random. Therefore, the most probable P is determined 
by the balance of the number of states and the weight factor, and the value of plaquette variable distributes around 
the most probable value for each point and each configuration. 

We first consider the case that there is no spatial correlation between the plaquette variables at each point and the 
volume is sufficiently large. In this case, the shape of the probability distribution as a function of the plaquette averaged 
over the space must be a Gaussian function. The central limit theorem tells us that the probability distribution of 
the average of the random numbers which have the same probability distribution is always of Gaussian type when 
the set of random numbers is large enough. We can apply this theorem in this case. Hence, 



where (P) is the expectation value of P and xp is the susceptibility, 

(P)= J P ^(P) dP, XP = GNsitciiP - {P)f) = GNsito j{P- {P)fw{P) dP. (9) 

We expect that w{P) is of Gaussian type also for more general interacting cases when the correlation length is much 
shorter than the size of the system. If we divide the space into domains which are larger than the correlation length 
and average the plaquette variables in these domains, the averaged plaquettes can be independent for each domain. 
When the number of domains is large, the distribution function as a function of the plaquette averaged over space 
must be a Gaussian function. 

However, we do expect that the probability distribution function is not of Gaussian type for the following two cases. 
One is, of course, the case that the correlation length is not small in comparison to the size of the system because the 
above-mentioned argument cannot be applied. The other case is that the most probable values of plaquette is not 
unique. For this case, the whole space is separated into domains having different states, and the plaquette variables 
in each domain distribute around one of the most probable values of plaquette. Although, on the surface separating 
these domains, the most probable plaquette value may not be realized, the effect from the wall becomes smaller as 
the volume increases, since the effect from the wall increases as a function of the area of the wall. Consequently, the 
existence of the domain wall does not affect the probability in the infinite volume limit. The probability distribution 
function should then be flat in the range between these most probable values of P because the spatial average of P 
depends on the size of these domains but the probability does not change in this range. However, in a finite volume, 
the effect from the domain wall cannot be neglected, hence the distribution function has two peaks when the number 
of most probable values for P is two. Clearly the two exceptions discussed here correspond to the case at a second 
order phase transition point and at a first order phase transition point, respectively. 

Here, it is convenient to introduce the effective potential defined by 

V{P) = -\nw{P). (10) 

As discussed above, the distribution function is normally written as 'w{P) ^ exp{— (6A'sitc/2xp)(^' — (-P))^}- When 
one considers a Taylor expansion of V{P) around the minimum (P), where the slope of the potential dV/dP is zero, 
the effective potential is dominated by the second order term in the region near the minimum, i.e. the potential is 
a quadratic function in the vicinity of (P), and the second derivative (curvature) of V{P) at (P) is related to the 
plaquette susceptibility with 

dP^ XP ■ ^ ^ 

A second order phase transition point is characterized by the slope and curvature of the effective potential. The 
slope dV/dP and curvature d^V/dP^ become zero simultaneously at the critical point. As given in Eq. ([5]), xp is an 
indicator of fluctuations and diverges at a second order phase transition point in the thermodynamic limit. When the 
susceptibility xp becomes large in the vicinity of a second order phase transition point, the effect from the second 
order term of V{P) becomes small in comparison to the higher order terms, and then the distribution function deviates 
from a Gaussian function. On the other hand, in the case of a first order phase transition point, more than one peak 
exist in the distribution function. This means that there are points which give dV/dP = more than once, and the 
curvature of V{P) may be negative around the mean value of P. 
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At the end of this section, we should also discuss the relation between the plaquette distribution function and the 
fourth order Binder cumulant, 

. <'^-'^'''>. (12, 
((P-iP))')' 

which often is used to identify the nature of a phase transition (25j . The value of the Binder cumulant at a second 
order critical point depends on the universality class. In the case of a first order phase transition, assuming the 
plaquette distribution is a double peaked function, the Binder cumulant is estimated as 



J{P - {P))'^w{P) dP _ 

(/(p- (p))2w(p) dpy ~ (A^ 



- TTTZ ,„x2 ~ 7X2X2 ~ y^-^) 



where the distance between two peaks is 2A and is wider than the width of each peak. On the other hand, when 
the distribution function can be modeled by a Gaussian function for a crossover transition or at a normal point, the 
Binder cumulant is given by 

B,-. v^!ip-{p)f^~'''-''~'^''^^^dp ( /Fdv^\ // ^^^^ 



(^./^ ^{P - {P)Ye-^iP-^P)ydPy 



In a region where a first order phase transition changes to a crossover, the Binder cumulant changes rapidly from one 
to three. We expect to find such a region for full QCD at high temperature and density. 

In addition, the method of Lee- Yang zeros has been used to identify the nature of the phase transition. The 
relation between the plaquette distribution function and the scaling analysis of the Lee- Yang zero has been discussed 
in Ref. [l^. The scaling behavior of the Lee- Yang zero can be also explained by the plaquette distribution function. 
Hence, the distribution function of the plaquette plays an important role in the investigation of the nature of a phase 
transition. 

III. LATTICE QCD AT FINITE DENSITY 

The most difficult problem for lattice studies at finite baryon density is that the Boltzmann weight is complex when 
the chemical potential is non-zero. In this case, the Monte-Carlo method is not applicable directly, since configurations 
cannot be generated with a complex probability. A popular approach to avoid this problem is the reweighting method. 
We perform simulations at /ig — 0, and incorporate the remaining part of the correct Boltzmann weight for finite 
/ig in the calculation of expectation values. Expectation values (O) at (f3,fj,q) are thus computed by a simulation at 
(/3o, 0) using the following identity, 

/Q^Nfiln dot A/(M,)-ln dct M(0)) g6(/3-/3o)A'=itoP\ 

_ J ' {1^0,0) /I r^ 

\^/{l3,t^g) /„JVf(lndetM(M,)-lndetM(0))p6(/3-/3o)Af.iteP\ ' ^ ' 

\ /(/3o,0) 

This is the basic formula of the reweighting method. However, because lndetM(/ig) is complex, the calculations of 
the numerator and denominator in Eq. (jlSp becomes in practice increasingly more difficult for larger iiq . We define the 
phase of the quark determinant 6 by the imaginary part of A'f In det M(/ig). If the typical value of 6 becomes larger 
than tt/2, the real part of e*^ {— cos9) changes its sign frequently. Eventually both the numerator and denominator 
of Eq. (fT5|) become smaller than their statistical errors and Eq. (fT5|) can no longer be evaluated. We call it the "sign 
problem" . The sign problem becomes more serious when the volume is large and the quark mass is small [23l [23 | . 

A. Reweighting method for finite fJ,q/T 

Let us discuss the reweighting method for finite fig using the plaquette distribution function. Originally, the 
reweighting method was proposed using the distribution function (his togr am) in Ref. (21.] . and applications to the 
finite density QCD in this style have been discussed in Ref. [lO, HI, I27L [28| . 

Here and hereafter, we restrict ourselves to discuss only the case when the quark matrix does not depend on P 
explicitly, e.g. the standard Wilson and staggered quark actions, the p4- improved staggered quark action etc., for 
simplicity. The partition function can be rewritten as 



Z(/3,/i,)= J R{P,^^g)w{P,|3) dP, 



(16) 
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where w{P, f3) is defined in Eq. ([6]) at /ig = and R{P, fiq) is the reweighting factor for finite iiq defined by 

J-pu 5{P'-P){detM{0))N' - ^ ^ 

This R{P, jiq) is independent of /3, and R{P, fiq) can be measured at any (3 using the following identity, 

where (• ■ ■ )(/3 =o) means the expectation value at /ig = 0. In this method, all simulations are performed at fiq — 
and the effect of finite fig is introduced though the operator det M{nq)/ detM(O) measured on the configurations 
generated by the simulations at /ig = 0. The expectation value of 0[P] is given by 



/ 0[P]R{P, iiq)wiP, f3) dP 
jR{P,fiq)w{P,(3)dP 



Moreover, the weight factor for non-zero /ig is R{P, iJ.q)w{P, (3), and thus the effective potential is defined by 

V{P, ^lq) = - ln[i?(P, fiq)w{P, /?)] = - lni?(P, /ig) + V{P, /?, 0). (20) 

The shape of the effective potential can then also be investigated at non-zero /ig once R{P, ^iq) is obtained. 

However, there are two problems to calculate R{P,iJ,q). The first problem is that the calculation of the quark 
determinant detM(/j,g) is very expensive. With present day computer resources, the exact calculation of detM(/ig) 
is difficult except on small lattices. The second problem is the "sign problem". Because lndetM(/ig) is complex, the 
calculations of the numerator of Eq. p8)) becomes in practice increasingly more difficult for larger /ig . If the complex 
phase factor of the quark determinant Re[e'^] changes its sign frequently, the expectation value of /ig) becomes 
smaller than its statistical error and the calculation of — In i?(P, Hq) in the effective potential becomes impossible. 



B. Taylor expansion in terms of /iq/T 



To avoid the first problem, we perform a Taylor expansion in terms of /ig around /ig 
coefficients, as proposed in Ref. 3]. We expand lndetM(/ig) in a Taylor series. 



In 



deiMipq) 
det M(0) 



OC ^ 



(9" (In det M) 



n\ L d{^iq/Ty 



and calculate the expansion 



(21) 



The Taylor expansion coefficients are rather easy to calculate by using the stochastic noise method. Although we 
must cut off this expansion at an appropriate order of /Xg, this approximation is valid at low density and can be 
systematically improved by increasing the number of the terms. 

Here, we discuss the effect of a truncation of the expansion. To estimate the range of fiq/T where the approximation 
is valid, an analysis of the radius of convergence is useful. The radius of convergence for pressure p{T, fj,q) is studied 
in Ref. [Q]. When one performs a Taylor expansion for p/T^ — lnZ/{VT^), 



p{T,iiq) p{T,0) 



2^4 



rp4 



(22) 



the radius of convergence can be defined by 



p = hm Pi, 



Pi 



Cj+2 



(23) 



for i = 2, 4, 6, • • • , where the odd terms are zero because the partition function is an even function of /ig. The Taylor 
expansion converges in the range of ^q/T < p when we consider all order of the expansion coefficients. This radius 
of convergence determine the lower bound of the critical point. This means conversely that the upper limit of the 
application range must be below the critical point if we estimate thermodynamic quantities using the Taylor expansion 
coefficients of the pressure or \nZ. 
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However, this problem may be avoidable when we consider a Taylor expansion of the reweighting factor i?(P, /ig) in 
Eq. Uni), since the weight factor itself should not be singular even at the critical point. Therefore, we expect that the 
application range is not limited by the critical point and evaluations beyond the critical point is possible. The same 
discussion of the radius of convergence is possible for the reweighting factor \nR{P,i.iq). We define the expansion 
coefRcients by 

oo 

\nR{P,^^,) = Y.''^{P){^)\ (24) 

1=1 

where the odd terms should be zero again. The radius of convergence is 



Pi 



li+2 



(25) 



When we neglect terms higher than 0{fj,^) in the calculation of IndetM, Eq. (PT|) . the application range can be 

estimate by Hq/T ^ pn ■ Because this approximation does not affect calculations of ri{P) for i < n, the truncation 
error is negligible when the contribution from higher order terms is smaller than that from the lower order terms. In 
the range where Hq/T < the (n + 2)*'' order term |r„+2(Mg/r)"''"^| is smaller than the n"^ order term |r„(/iq/r)"|. 
Hence, the truncation error must be small in this range. 

Before discussing the radius of convergence for In R{P,fj,q) using the data obtained by Monte Carlo simulations, 
we estimate the application range at large P and small P, corresponding to large temperature and small temper- 
ature, respectively. In the free quark gas limit, where P is maximum, the quark determinant is expected to be 
{lndetM){Nt/Nsf = {Ttt^/GO) + (l/2)(/i,/T)2 + (l/47r2)(^q/T)4 in the continuum limit ^. Because is infinity, 
the convergence of the Taylor expansion seems to be good for large P. 

On the other hand, in the study of the equation of state J^, |20], the numerical results of the derivatives of pressure 
with respect to fiq/T at low temperature have been found to reproduce the prediction from the hadron resonance gas 
model very well. Because small plaquette values are generated in the low temperature simulation, this model may 
give a suggestion of the application range for small P. The quark chemical potential dependence of pressure in the 
hadron resonance gas model is discussed in Ref. It is suggested that 

Pipq) p(0) ( 3pq 



rp4 rp4 1 ji 

and the radius of convergence for pressure is given by 



(X cosh {-^ , (26) 



This Pi increases as i increases, and the convergence radius p is infinity. Although we expect that the radius of 
convergence for InR is larger than that for the pressure, we try to estimate the application range from this pi. When 
we neglect terms higher than 0{iJ,q), as it is done in this study, the application range is suggested to he pq/T < pQ w 2.5. 
This implies that the error that arises from the approximation up to 0{pq) may be sizeable for Pq/T 2.5, and more 
careful arguments are required when we calculate the reweighting factor R{P, Pq) for Pq/T 2.5. We will discuss this 
application range in Sec. IIII Fl again. The results of In R{P,pq) obtained by the calculations up to 0{pq) and 0{pq) 
will be compared, and we will confirm that the truncation error is still small even at Pq/T 2.5. 



C. Avoidance of the sign problem at finite density 



We discuss here how to avoid the sign problem in our reweighting approach. In the framework of the Taylor expan- 
sion, we can easily separate \nAet M{pq) into real and imaginary parts because the even derivatives of IndetM(^q) 
are real and the odd derivatives are purely imaginary [3]. The absolute values of the quark determinant and the 
complex phases 6 are thus given by 

N, In I det A/| = ]V,Ite |l.(del M)| = N, f ^Re^3=i|i^ (^l) , (28) 
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e=(A?^4)Im[ln(det M)] 



FIG. 1: The histogram of the complex phase for /ig/T = 1.0 and 2.0 at /3 = 3.65 {T/Tc = 1.00) on a 16^ x 4 lattice. The dashed 
lines are the fit results by Gaussian functions. 



where one must replace N{ in these equations to 7Vf/4 when one uses a staggered type quark action. Here, it is worth 
noting that 9 corresponds to the complex phase of the quark determinant however by definition this quantity is not 
restricted to the range from — tt to tt because there is no reason that the imaginary part of Indet M in Eq. ((29)) must 
be in the finite range. In fact, this quantity becomes larger as the volume increases. 

We show histograms of 6 at the pseudo-critical temperature (/3 = 3.65) for fiq/T =1.0 and 2.0 in Fig. [I] where 9 is 
calculated using the data of the Taylor expansion coefficients up to 0{fi^) obtained with two- flavors of p4-improved 
staggered quarks in Ref. 0. These histograms seem to be almost Gaussian functions. We fit these data by Gaussian 
functions, ^ exp{—x9'^), where the overall factor and x are the fit parameters. The dashed lines in Fig. [T] are the ffi 
results. It is found that the histogram of 9 is well-represented by a Gaussian function. 

Similar to the discussion of the Gaussian distribution function for the plaquette in Sec. [ill we may argue that 
the histogram of is a Gaussian function. Because there is no critical point in two-flavor QCD with finite quark 
mass at fiq — 0, the spatial correlation length between the quark fields is not expected to be long. The Taylor 
expansion coefficients in Eq. (P^)) are given by combinations of traces of products of 9"A//9(/ig/T)" and (see 
Appendix of Ref. [3]). Therefore, the expansion coefficients are obtained by the sum of the diagonal elements of 
such matrices. When the correlation among the diagonal elements is small and the volume is sufficiently large, the 
distribution functions of the expansion coefficients and 9 should be of Gaussian type due to the central limit theorem. 
For example, the diagonal elements of the ffi'st coefficient, Im[(9(lndet M)/9(/ig/T)] — Im[tr[M~^(9M/9(/ig/T))]], 
is the imaginary part of the local number density operator at /i^ = 0. If the spatial density correlation is not very 
strong, the Gaussian distribution is expected. Figure [1] is consistent with this argument. 

We note that, once we assume a Gaussian distribution for 9, the problem of complex weights can be avoided. 
A variety of distribution functions with respect to various quantities are discussed in the density of state method 
[20l . [2ll . [26l . [27I . [28j . We introduce the probability distribution iD as a function of the plaquette P, the absolute value 
of [detM(/iJ/detM(0)]^' = F and the complex phase ee Im[lni^(^g)], 

w{P\ \F\', 9') = J VU5{P' - P)5{\F\' - \F\)5{e' - 9)(dei Af(0))^feS'^^=""^. (30) 

The distribution function itself is defined as an expectation value at jiq = 0, i.e. w{P', \F\',9') cx {S{P' — P)6{\F\' — 
\F\)6{9' — d)){T,fiq=o)j however \F\ and 9 are functions of /ig/T obtained by the Taylor expansion at /ig = 0. The 
expectation value of 0[P, \F\, 9] a,t fig = is given by 

- ^(^17^ /^^/ '^l^l / (^[PAnd] ^{PAF19). (31) 

Since the partition function is real even at non-zero density, the distribution function has the symmetry under the 
change from 9 to —9. Therefore, the distribution function is a function of 9^, e.g., w{9) ~ exp[— (026'^ -|- 04^'* -|- Og^^ -|- 
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■•■)]. Moreover, as we discussed, when the system size is sufSciently large in comparison to the correlation length, 
the distribution function should be well-approximated by a Gaussian function: 



w{P, \F19) « J |F|) exp [-a2(P, \F\)e'] . (32) 



We assume this distribution function in terms of 6 when P and \F\ are fixed. 
The coefficient a2{P,\F\) is given by 



1 



2a2{P',\F\') 



J de 6*2 w{P', \F\', B) I J de w{P', \F\', 0) 
{eH{P'-P)6{\F\'^\F\))^^^^^^^^ 
{5{P'-P)5m'-\F\))^T,,,l,) ' 



(33) 



using ^/a^ J 0^ exp{-a2e^)de = l/(2a2). 

When the volume is sufficiently large, this assumption will be valid except at a critical point. For the case of 
two-flavor QCD at finite quark mass, this assumption should be valid because there is no critical point for Hg = 
except in the chiral limit, and is suggested by Fig.[Tl though the values of P and |F| are not fixed in the calculation 
of Fig. [TJ Then, the integration over 6 can be carried out easily and we obtain the numerator of Eq. p8|) for the 
calculation of R{P, fig), 

{F{^^MP'-P))^T,,,=o) ^j'^P j ^1^1 / de^w'{P^F\)e-^^^o\^'\F\5{P' -P) 

= ^ JdP J d\F\ w\P,\F\)e-^^^^'''^\F\S{P' - P) 

= ^ y"l?[/e-i/(4-^(^^l^l»|F(Ai,)|(5(P'-P)(detA/(0))^'e-^» 

= (e-V(4a.(P.|i^|))|^(^^)|^(p._p)\^ (34) 



Since 9 is roughly proportional to the size of the quark matrix M, the value of l/a2 becomes larger as the volume 
increases. Therefore, the phase factor in R{P, /ia) decreases exponentially as a function of the volume. However, 
the most important point in this approach is that the operator in Eq. p4p is always real and positive for each 
configuration in this framework, hence the expectation value of R{P, /ig) is always larger than its statistical error, i.e. 
the contribution In R{P, fig) to the effective potential V{P, P, fig) is always well-defined. Therefore, the sign problem 
is completely avoided if we can assume the Gaussian distribution of 6. 

We calculate 6 using the stochastic noise method. Then, the value of 6 contains an error due to the finite number of 
noise vectors (iVn). As discussed in Ref. [3, 2^, a careful treatment is required to reduce this error for the calculation 
of yJJW), i.e. width of the distribution of 9. Since the noise sets for the calculation of the two 9 in the product must 
be independent, we subtract the contributions from using the same noise vector for each factor. By using this method, 
we can make the TVn-dependence of {9"^) much smaller than that by the naive calculation from rather small iVn, 
hence it may be closer to the = oo limit. We took iVn = 50 or 100 in this calculation, so that the iVn-dependence 
is negligible. On the other hand, as increases, the result of obtained by the naive calculation without the 

subtraction becomes smaller and approaches the result with the subtraction. For the case at /3 = 3.65, fiq/T = 2.0 
with Nn = 100 in Fig. [1] the di fferen ce between them is about 13%. Since the width of the distribution function 
shown in Fig. [1] corresponds to -y/ (9^) without the subtraction, the width in Fig. [T] is slightly larger than that in the 
A^n = oo limit. 

For more quantitative arguments of the Gaussian distribution function, we also compute the fourth order Binder 
cumulant of the complex phase for fiq/T — 1.0 and 2.0, using the data obtained in a simulation of two-fiavor QCD 
with p4- improved staggered quarks, Ref. 7]. The Binder cumulant is defined by 

Bi.^^ (35, 

As discussed in Sec. |lTl this quantity is a good indicator to check whether the distribution is of Gaussian or not. To 
confirm the validity of the assumption, Eq. ([5^ . we should compute as a function of P and However, because 
the width of the plaquette distribution function w(P, /3) for each (3 is narrow in our simulation (see the results in 
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FIG. 2: The fourth order Binder cumulant of the complex phase for ^iq/T = 1.0 (circle) and 2.0 (square). The filled symbols 
are the results obtained when the contributions from using the same noise vector are subtracted in the products of 9. The open 
symbols are the results without the subtraction. The dashed line is the value of Gaussian distribution. 



Sec. IIII D|) . we calculate for each j3 without separating the configurations in terms of P. The circle and square 
symbols in Fig. [2] are the results for Hq/T = 1.0 and 2.0, respectively. We use the stochastic noise method for the 
calculation of the products of 9. The results plotted by filled symbols are obtained when the contributions from using 
the same noise vector are subtracted. The open symbols are the results without the subtraction. Because the complex 
phase vanishes in the large (3 limit, {9"^) becomes smaller as (3 increases. We omitted results having large statistical 
errors due to the small {9"^) at large (3. We find from this figure that the results of i?f are almost consistent with three. 
As we discussed in Sec. if the distribution is of Gaussian, the Binder cumulant is three. Hence, this figure suggests 
the Gaussian distribution. The results around (3 = 3.66 are slightly larger than three, but the difference would be 
within the systematic error due to finite statistics because usually the Binder cumulant becomes smaller than three 
when the correlation length is long. 

To estimate the effect when the distribution is slightly different from Gaussian, we consider a distribution func- 
tion with small 04, i.e. 'w{9) ~ exp[— 02^^ — 046''*]. In this case, the phase factor changes from exp[— l/(4a2)] to 
exp[— 1/(402) + 804/(402) — 04/(160!) + •••], and also the expectation value of 9"^ for fixed P and |F| becomes 
(6*^) = 1/(202) — 804/(203) + • • • . Since the term of 804/(402) is absorbed into (0^)/2, the leading contribution from 
04 in the phase factor is exp[— 04/ (1602)]. Because 1/02 ~ 0(/x^), this effect becomes larger as increases. Therefore, 
for the case of 04 ^ (Bf 7^ 3), the estimation of the range of jiq in which the non-Gaussian contribution is small 
may be important as well as the application range of the Taylor expansion in /i^ discussed in the previous section. 



D. Reweighting method for (3 direction 



We consider the reweighting method for the /? direction at ^iq = 0. This is the case of R{P,0) — 1. Using the 
plaquette distribution function (plaquette histogram) w{P,Po) at the simulation point /3o, the expectation value of an 
operator given by the plaquette is evaluated by 

where we discuss only the case when the quark matrix does not depend on (3 explicitly for simplicity, otherwise 
equation (|36p is no longer correct. 

From Eq. ([55]) . under the parameter change from (3o to (3, the weight w{P, (3) becomes 

w{P, 13) = eS('3-''°)^='*=^u;(P, Po). (37) 
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FIG. 3: The plaquette histogram and the effective potential at = as a function of the plaquette for the two-flavor 
p4-improved staggered action obtained in Ref. 0- 
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FIG. 4: The curvature of the effective potential at /i, — 0. 



If we rewrite e ^^''^''*''^w{P, (3o) = w{P), we obtain Eq. ([T]) from Eq. p6|) . The effective potential becomes 

V{P, 13)^^ In u;(P, /3) = V{P, M - 6(/5 - [io)N,it.P- (38) 

When (3 is increased, the slope of V{P) becomes smaller, whereas the curvature of V{P) does not change. This 
implies that the curvature of V{P) is independent of j3. For the case of cPV/dP'^ > 0, the value of P which gives the 
minimum of V{P) becomes larger as (3 increases. 

Here, we want to explain the /3 dependence of the effective potential using the data from Ref. .2]. The configurations 
were generated with Symanzik-improved gauge and two- flavor p4-improved staggered fermion actions. Because the 
improved gauge action was used in Ref. f7'| , the definition of P is 



P = -- 



Sq 



(39) 



where W^^"^ is the I x J Wilson loop for each point and each direction. The maximum of this P is 1.5. 

The probability distribution function w{P), i.e. the histogram of P, and the effective potential V{P) are given in 
Fig. O These are measured at sixteen simulation points from (3 = 3.52 to 4.00 for the bare quark mass ma = 0.1. The 
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FIG. 5: The reweighting factor R{P,fiq) for fiq/T = 0.5 - 2.5 obtained by the Taylor expansion up to 0(/i^). The dashed lines 
are the cases when the effect of the complex phase is omitted, R{P,^q). 



corresponding temperature normalized by the pseudo-critical temperature is in the range oiT/T,. = 0.76 to 1.98, and 
the pseudo-critical point [T/T^ = 1) is about (i = 3.65. We show the values of (i and T/T(. above these figures. The 
ratio of pseudo-scalar and vector meson masses is mps/my fa 0.7 a.t /3 — 3.65. The lattice size Nsn^ is 16^ x 4. The 
number of configurations is 1000 - 4000 for each p. Further details on the simulation parameters are given in Ref. 0]- 
To obtain w{P) and V{P), we grouped the configurations by the value of P into blocks and counted the number of 
configurations in these blocks, and the potential V{P) is normalized by the minimum value for each temperature. 

Because the transition from the hadron phase to the quark-gluon plasma phase is a crossover transition for two- 
flavor QCD with finite quark mass, the distribution function is always of Gaussian type, i.e. the effective potential 
is always a quadratic function. The value of the plaquette at the potential minimum increases as /3 increases in 
accordance with the above argument. 

Figure [3] shows the curvature of the effective potential at /i^ ~ 0, <PV / dP^ {P) = — (P {In w)/dP'^, as a function of P. 
We estimate this quantity from the relation between the plaquette susceptibility xp and the curvature of the potential 
at = 0, Eq. pT|) . Here, it should be emphasized again that the slope of the potential changes as Eq. ((38)) when (3 
is changed, but the curvature of the potential never changes. This means that the curvature is independent of /3 and 
is determined by the measure DU and the quark determinant det M. As we discussed in Sec. |lTl the curvature of the 
effective potential V{P) at P for the potential minimum is important to categorize the nature of phase transition, e.g. 
the curvature must be zero at a second order phase transition point. The property of the curvature being independent 
of P will make our analysis simpler in the next section. 



E. Numerical calculations of the reweighting factor 



We calculate the probability distribution function at non-zero fiq using the data of the Taylor expansion coeffi- 
cients up to 0{iig) computed in Ref. 0] with the p4-improved staggered quark action. Since the simulations are 
performed in the region where no critical points exist, the assumption of the Gaussian function is valid. The coef- 
ficient a2{P, \F\) in the distribution function of 9 is measured using Eq. ([551) . However, because the values P and 
|F| = I det M (fiq) / det M {0)\^^ on each configuration are strongly correlated [2J], we may assume that |F| is approx- 
imately given as a function of P for each configuration so that 02 (P, \F\) is given by a function of P only. In this 
approximation, the contribution from the complex phase in R(P' , fiq) can be simplified, 

Although the correlation between \F\ and a2 is neglected in this equation, the main contribution to the variation 
of R(P, fiq) comes from \F\, and the contribution from the phase factor is not very large, as we will see in Fig. [S] 
Therefore, the correlation of these two factors is negligible in the following argument. For the calculation of R{P, fiq), 



12 



8000 



6000 



4000 



2000 



I I 

d{lnR) 


1 1 




III 


dP 














n /r=2.5 

LiV=2.0 


J 






LI /r=i.5 

n/r=i.o 








^''/r=o.5 










1,1, 


1 , 1 




1,1,1, 



0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 



FIG. 6: The slope of \n R{P, fiq) as functions of the plaquette. 



we use the delta function approximated by a Gaussian function, d{x) « l/(A^/7r) cxp[— (x/A)^], where A = 0.0025 is 
adopted. 

Because R{P,fiq) is independent of /3, we mix all data obtained at different /3. This mixture can be justified by 
extending Eq. (0 for multi-/3, e.g. R{P' , ^iq) = [Ni{5{P' - P)F) N2{5{P' - P)F) p,]/[N^{5{P' - P)) N2{5{P' - 
P))^^] for the data at (3i and [32 with the number of configurations iVi and N2. The results for In i?(P, ^g) are shown 
by solid lines in Fig. [5] for /i^/T — 0.5, 1.0, 1.5, 2.0 and 2.5. We find a rapid change in Ini? around P ~ 0.83, and the 
variation becomes larger as [lajP increases. 

The dashed lines in Fig.[5]are the results that we obtained when the effect of the complex phase, i.e. exp[— l/(4a2)], 
is omitted. We define this quantity as 



i?(P',/i,) 



(|p(m,)|^(p'-p)) 



('5(^''-P))(T..,=, 



(41) 



0) 



We discuss in Sec. IIII Gl that these dashed lines correspond to the reweighting factor with non-zero isospin chemical 
potential /x/ and zero quark chemical potential //g, i.e. R{P, fiq) — P(P, ///). The variation of InP in terms of P 
becomes milder when the effect of the complex phase is omitted. 

The effective potential V{P, f3, Hq) is obtained from Eq. ([^0]) substituting the data in Fig. [3] and Fig. (5] To study 
the existence of a second order phase transition, the curvature of the potential is important. The minimum of the 
potential can be changed by shifting /? but the curvature can be controlled only by \n R{P, fiq). The result of the 
curvature at fiq = 0, — (P (in w)/dP'^, as a function of P is shown in Fig. U) Because — (P (In ■w)/dP'^ is positive, a 
region where (P {In R) / dP^ > is necessary for the existence of a critical point. The curvature of InP is positive for 
P < 0.83. 

In order to analyze the sign of d^V/dP^{P, fiq), we fitted the data around P by a quadratic function, lnP(P', fiq) — 
xo+xi{P' ^P)+X2{P' — P)^, where xq, xi and X2 are the fit parameters, and calculate the first and second derivatives 
oflnP(P,^q) at each P. The result of the slope, d([nR)/dP{P, fiq) = a; i, is shown in Fig. [6] for each ^q/T. We adopt 
the result obtained by fitting in the range between P — 0.015 and P + 0.015 for each P as the final result. In the region 
around P ~ 0.83, d(lnR)/dP becomes larger as fiq/T increases and \nR{P,fiq) changes sharply in this region. The 
result of the curvature, d^ (In R) / dP^ {P, fiq) = 2x2, is plotted as solid line in Fig. [71 The magnitude of the curvature of 
InP also becomes larger as fiq/T increases. The dashed line in Fig.[7]is the data of —cZ^ (In u')/dP^(P) in Fig. ID This 
figure indicates that the maximum value of d^ (in R) / dP^ {P, fiq) at P = 0.80 becomes larger than — d'^ (In w) /dP"^ for 
fiq/T > 2.5. This suggests that the curvature of the effective potential, d'^V/dP^ = - d"^ (In w)/dP^ - d"^ (In R) / dP^ , 
vanishes at fiq/T ~ 2.5 and a region of P where the curvature is negative appears for large fiq/T. 

Next, we estimate the value of /3 which gives the potential minimum at P = 0.8 for fiq/T — 2.5 by solving the 
equation: 



(42) 
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FIG. 8: The radius of convergence, p2, pi, for the Taylor expansion of R{P, pq) 



This equation can be solved without changing fJ^q/T, and tells us the location of the critical point in the iJ-q/T) 
plane. Since a simulation with (3 w 3.56 gives d{lnw)/dP = at P w 0.8, we adopt (3o = 3.56. Substituting 
d{\nR)/dP « 4000 at {P,iJ,q/T) = (0.8,2.5) in Fig. [6] and TVsitc = 16^ x 4, we obtain P « 3.52. This P corresponds to 
T/Tc = 0.76, where Tc is the pseudo-critical temperature at fiq = 0. Therefore, it is found that the potential is flat 
up to second order in P around P = 0.80 with {T/Tc, ^iq/T) k, (0.76, 2.5), suggesting the existence of a critical point 
around this value. 

Further studies are, of course, needed for the precise determination of the critical point in the {T,fiq) plane, 
increasing the number of terms in the Taylor expansion of Indet M and decreasing the quark mass in the simulation. 
The quark mass is still heavier than the physical quark mass. However, the arguments given above indicate the 
existence of a first order phase transition line at large ^iq/T because the magnitude of the curvature of R{P, /iq) 
increases monotonically and eventually the curvature of the potential becomes negative at large fJ-q/T, corresponding 
to a double-well potential of a first order phase transition. 
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FIG. 9: The reweighting factor R{P,fiq) computed by the Taylor expansion up to 0(/i^) (sohd hues) and 0{fig) (dashed lines). 



F. Application range of this analysis 

Next, we discuss the reliability of our analysis in view of the truncation of the Taylor expansion used here. Because 
the dominant contribution in Ini? is given by the reweighting factor without the phase effect, Ini?, we consider the 
radius of convergence for In R. The expansion is defined by 

oo 

lni?(P,^,) = 

n=l 

r2 = {d2)p, fi = {di)p + i {{dl)p - {d2)l) , 

fe = {d6)p + {d2di)p - {d2)p{di)p + i {{dl)p - 3{d2)p{dl)p + 2{d2)%) , (44) 

where d„ = (A^f/n!)a"(lndet Af)/(9(/^g/T)", {■■■)p' = {■ ■■d{P' - P))/{S{P' - P)), and the odd terms are zero. The 
radius of convergence is obtained by analyzing the asymptotic behavior of pn — VrnV^n+al for n = 2, 4, 6, • • • , oo. 

In this analysis, we calculated Indet A/ using the data of dn up to 0(/iq). This approximation does not affect the 
calculations of f2,r4 and fg, but there is a missing term, i.e. {ds)p, in the calculation of fg. If the 8**^ order term of 
Ini? is larger than the 6'^ order term, the effect of the truncation may be sizeable. Because \ fQ{p,q/T)^\ > |f8(/Xq/r)^| 
for fiq/T < pq, the application range for our current analysis should be Pq/T < pQ. We calculate p2 and p^. These 
results are shown in Fig. [H] The dashed line is p2 in the free gas limit, and p4 is infinity in the free gas limit. The 
results of r2 and are positive for all P we investigated, but fg changes its sign at P = 0.84. fg is negative for 
P > 0.84. We find that p4 (square) is larger than p2 (circle), and the values of p2 and p4 are larger than the hadron 
resonance gas model values, p2 ~ 1-15 and p4 w 1.83. For our analysis, where we omitted the calculation of dn higher 
than the 6**^ order in Ini?, the application range given by pg would be larger than the hadron resonance gas model 
prediction, pQ « 2.49, and the parameter range we investigated thus seems to be within the application range. 

We moreover estimate the effect from higher order terms in the Taylor expansion by changing the number of terms 
in the Taylor expansion. Figure M shows the difference between the results up to 0{pq) and 0{p,q) for Pq/T — 
0.5, 1.0, 1.5,2.0 and 2.5. The dashed lines are the same as the solid lines in Fig. [S] and the solid lines are the results 
obtained when the highest order term and the next highest order term, d^ (in M) / d(pq /T)^ and d^ (in M) / d{pq /T)^ , 
are omitted in Eq. (pij) . It is found from this figure that the difference becomes visible at Pq/T ~ 2.5, but the 
truncation error of the Taylor expansion does not affect the qualitative argument of the effective potential at finite 
density in the range we have discussed. For more quantitative investigation of the critical point in the (T, pq) plane, 
more accurate calculations including higher order terms in the Taylor expansion of pq may be important. 
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G. QCD with an isospin chemical potential 



Finally, it is worth discussing the difference between QCD with a quark (baryon) chemical potential and an isospin 
chemical potential. The isospin chemical potential is defined by fij = (/i^, — /Zd)/2, where and fj-d are the chemical 
potential for u and d quarks, respectively. For the case with non-zero isospin and zero quark chemical potentials, 
/^g — ifJ'u + M(i)/2 =0, i.e. /i„ = —^id — A*/, the quark determinant is real and positive because 



det M(^„)detM(/Xd) = det M (pi) det M{~^ii) = det M(^/)(det Af(/^/))* = |detM(/i/)|^ 



(45) 



where we used an identity at finite fXq, j5M{iJ.q)-f5 = M{—fj,qy. Therefore, Monte-Carlo simulations are possible for 
this case [s^, and the simulations with the isospin chemical potential have been performed in Ref. [isl. [T9l.l28l. Isil. Isl] . 
It may be important toward the understanding of QCD at finite density to consider the difference between the phase 
diagram with non-zero baryon chemical and that with non-zero isospin chemical potential, 

The reweighting factor R, i.e. the dashed line in Fig.[5l corresponds to the reweighting factor of the isospin chemical 
potential R{P, /i/) for each /i//T because the quark determinant is | det It is found from Fig.[5]that the slope 

and the curvature of Ini? around P ~ 0.82 for the isospin chemical potential are smaller than those for the quark 
chemical potential. This means that the value of jJLi/T where the second order phase transition appears by canceling 
the curvatures of lnw(P, /3) and \nR{P,iJ,i) is larger than the critical point of Hq/T. It is suggested in Ref. [l3| that 
there is no first order phase transition region in the low density regime of QCD with non-zero fii/T. Although more 
quantitative estimations of the reweighting factor are needed to confirm the existence of the first order transition line, 
our argument may be related to their result. 

Furthermore, in the case of the approximation up to 0{^^ j), R{P,iiq) and R{P,^j) have a close relation to the 
quark number susceptibility Xg and isospin susceptibility xi at Mg,/ = 0- Using the equations (|28p . (|29p and 



/ fl d^(\ndetM) //i„\2 
lni?(P,Mg)«ln(exp|-MRe^^-^(^^ 



Nf 



92 (In dot Af) 
5(Mg/T)2 



Nf 



9(lndet A/) 



in this approximation, and when the effect from 9 is omitted, we find 



1// 9(lndetM)Mg 



, . X , ^ 1 / d^(\ndetM) 

lnfl(P.M,) = l„«(P,„)-5(A',-i^^ 




(46) 



(47) 



where (• • • )p' = (• • • S{P' — P)) / {S{P' — P)). These are related to Xq/T'^ and Xi/T'^ as functions of /? (temperature) 
by the following equations 



J^{T,fIq,I^O) 



Nll_ 

mz 



Xi_ 

T2 



m/^g./ = 0) 



a2(lndetAf) 
92(lndetAf)\ 



/ 9(lndetA/) ' 
w{P,(3) dP. 



w{P,(3) dP, 



(48) 
(49) 



From these equations, the similarity between the figures for R(P, fj,qj) and those of the quark number and isospin 
susceptibilities can be easily understood in the regime where the Taylor expansion is valid. As shown in Fig. [31 w{P, /3) 
is a Gaussian function having a sharp peak. Therefore, Fig. [5] is quite similar to Fig. 1 in Ref. [7] if we replace the 
horizontal axis P by T/Tc{P). As we have discussed, the positive curvature in the P dependence of lni?(P, /u,) is 
required for the appearance of the critical endpoint. It is found that the positive curvature is related closely to the 
rapid increase of the quark number susceptibility near the pseudo-critical temperature at /ig = 0. 

Here, it should be noted that Xi/T"^ is always larger than Xq/T"^ at /ig = because 9(lndet Af )/9(/ig/T) is purely 
imaginary, and thus (9(lndet Af)/9(/Xq/T))2 is negative. Moreover, both susceptibilities approach the same value in 
the high temperature limit. Hence, the variation of xz/T^ around the transition point would be milder than that of 
Xq/T"^, corresponding to the behavior of R{P, fiq) and R{P,fii). This may explain the difference between the phase 
diagrams with finite /ig and finite (Xj. Furthermore, in the framework of the hadron resonance gas model at low 
temperature, the isospin susceptibility correspond to fluctuations of pions, and the pion mass is more sensitive to the 
quark mass than baryon masses. Therefore, when the quark mass is decreased, the pion mass becomes smaller and 
the fluctuation becomes larger at low temperature. This suggests the change of Xj/T^ around Tc may be milder at 
small quark mass, i.e. the difference between lni?(P, /Zg) and In R{P, fij) becomes large at small quark mass. 
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For more precise arguments on the phase structure, more accurate evaluations of R{P, ^iq) and /i/) are required 
increasing the number of terms in the Taylor expansion of In det M . However, the qualitative property that the critical 
value of ^lq/T in the {T,iiq) plane is smaller than the critical /ij/T in the {T,^j) plane can be understood by the 
well-known properties of the quark (baryon) number and isospin susceptibilities combined with the argument of the 
effective potential. 

IV. CONCLUSIONS 

We have discussed the phase structure of lattice QCD at non-zero density. The probability distribution as a function 
of the plaquette was estimated at non-zero temperature and chemical potential using the data obtained with two- 
flavors of p4-improved staggered quarks in Ref. Q. In this analysis, we have adopted two approximations. One is 
that we estimate In det M from the data of a Taylor expansion up to 0(/Lt^). Terms of higher order than are 
omitted. We have estimated the range where this approximation is valid and studied in the reliability range. The 
second approximation is an assumption on the probability distribution for the complex phase. We have assumed the 
distribution function to be a Gaussian function. This assumption will be valid for sufficiently large volume and we 
have checked that the distribution is well- approximated by a Gaussian function for the data used in this analysis. 

In spite of the use of these approximations, it is found that the shape of the effective potential which is of Gaussian 
type at Hq = changes to a double-well type at large iJ-q/T. This property is related closely to a well-known behavior of 
the quark number susceptibility at iiq — 0, i.e. the rapid increase near the phase transition point. For the quantitative 
estimation of the endpoint of the first order phase transition, further investigation must be needed. However, this 
argument strongly suggests the existence of the first order phase transition line in the (T, fj,q) plane. 

We also discussed the difference between QCD with a quark chemical potential and QCD with an isospin chemical 
potential, and found that the critical value of the quark chemical potential seems to be smaller than that of the isospin 
chemical potential. 
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